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Abstract. This paper gives a new representation of Pickands' constants, which arise in the 
study of extremes for a variety of Gaussian processes. Using this representation, we resolve 
the long-standing problem of devising a reliable algorithm for estimating these constants. A 
detailed error analysis illustrates the strength of our approach. 



1. Introduction 

Gaussian processes and fields have emerged as a versatile yet relatively tractable class of 
models for random phenomena. Gaussian processes have been applied fruitfully to risk theory, 
statistics, machine learning, and biology, while Gaussian fields have been applied to neu- 
roimaging, astrophysics, oceanography, as well as to other fields. Extremes and level sets are 
particularly important in these applications [7j- New applications and theoretical developments 
continually revive the interest in Gaussian processes, see for instance |27j . 

Although the understanding of Gaussian processes and fields has advanced steadily over the 
past decades, a variety of results related to extremes (tail asymptotics, extreme value theorems, 
laws of iterated logarithm) are only 'explicit' up to certain constants. These constants are 
referred to as Pickands' constants after their discoverer [29]. It is believed that these constants 
may never be calculated [lj. 

These constants have remained so elusive that devising an estimation algorithm with certain 
performance guarantees has remained outside the scope of current methodology [18]. The 
current paper resolves this open problem for the classical Pickands' constants. Our main tool is 
a new representation for Pickands' constant, which expresses the constant as the expected value 
of a random variable with low variance and therefore it is suitable for simulation. Our approach 
also gives rise to a number of new questions, which could lead to further improvement of our 
simulation algorithm or its underlying theoretical foundation. We expect that our methodology 
carries through for all of Pickands' constants, not only for the classical ones discussed here. 

Several different representations of Pickands' constants are known, typically arising from 
various methodologies for studying extremes of Gaussian processes. Hiisler [23] uses triangular 
arrays to interpret Pickands' constant as a clustering index. Albin and Choi [3j have recently 
rediscovered Hiisler's representation. For sufficiently smooth Gaussian processes, various level- 
crossing tools can be exploited [61 [25]. Yet another representation is found when a sojourn 
approach is taken [9J. Aldous [4J explains various connections heuristically and also gives intu- 
ition behind other fundamental results in extreme-value theory. We also mention Leadbetter 
et al. |26t Ch. 12], who use methods different from those of Pickands but arrive at the same 
representation. 

The approach advocated in the current paper is inspired by a method which has been applied 
successfully in various statistical settings, see |33] and references therein. This method relies 
on a certain change-of- measure argument, which results in asymptotic expressions with a term 
of the form K(M/S), where M and S are supremum-type and sum-type (or integral-type) 
functionals, respectively. This methodology can also be applied directly to study extremes of 
Gaussian processes, in which case it yields a new method for establishing tail asymptotics. 
This will be pursued elsewhere. 
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Throughout this paper, we let B = {Bt : t G M} be a standard fractional Brownian motion 
with Hurst index a/2 £ (0,1], i.e., a centered Gaussian process for which 

1 



Cov(B s ,B t 



[\s\ a + \t\ a - \t 



Note that has stationary increments and variance function V&r(Bt) = \t\ a . The process {Zt} 
defined through Z t = V2B t — \t\ a plays a key role in this paper. This stochastic process plays a 
fundamental role in the stochastic calculus for fractional Brownian motion [8]. The 'classical' 
definition of Pickands' constant % a is 



(1.1) 



U a = lim ^E 

T-s-oo T 



sup 

te[o,T] 



Current understanding of % a and related constants is quite limited. It is known that Hi = 1 
and that H2 = 1/ PiS SO] , and that 7i a is continuous as a function of a [H] . Most existing 
work focuses on obtaining sharp bounds for these constants [U H5j H7J HHJ E21 [22]. Previous 
work on estimating Pickands' constant through simulation has yielded contradictory results 

PI EE]. 

The next theorem forms the basis for our approach to estimate H a - Note that the the- 
orem expresses T~L a in the form E(M/S). A different but related representation is given in 
Proposition [2] below, and we give yet another representation in Proposition [4] 



Theorem 1. We have 



E 



sup te R e 



z t 



tdt 



The representation E(M/S I ) is well-suited for estimating Pickands' constant by simulation. 
Although both M and S are finite random variables with infinite mean, we provide theoretical 
evidence that their ratio has low variance and our empirical results show that this representa- 
tion is suitable for simulation. 

This paper is organized as follows. Section [2] establishes two results which together yield 
Theorem [TJ In Section |3} we state an auxiliary result that plays a key role in several of the 
proofs in this paper. Section [4] gives an error analysis when K(M/S) is approximated by a 
related quantity that can be simulated on a computer. In Section [5j we carry out simulation 
experiments to estimate Pickands' constant. Some proofs are deferred to Appendix A, and a 
table with our simulation results is included as Appendix B. 



2. Representations 

This section is devoted to connections between Pickands' classical representation and our 
new representation, thus establishing Theorem [TJ We also informally argue why our new 
representation is superior from the point of view of estimation. This is explored further in the 
next section. 

The following well-known change-of- measure lemma forms the basis for our results. 

Lemma 1. Fix i £ 1, and set = {\/2B s — \s — t\ 2H : s £ M.}. For an arbitrary measurable 
functional F on IR K , we have 

Ee Zt F(Z) = EF(\t\ 2H + Z (t) ). 

When the functional F is moreover translation-invariant (invariant under addition of a con- 
stant function) , we have 

Ee Zt F(Z) = EF(9 t Z), 
where the shift Ot is defined through (0tZ) s = Z s —t- 

Proof. Set Q(A) = E[e Zt lA], and write for the expectation operator with respect to Q. 
Select an integer k and s\ < S2 < • • • < Sk- We show that (Z Sl , . . . , Z Sk ) under Q has the same 



ON ASYMPTOTIC CONSTANTS IN THE THEORY OF EXTREMES FOR GAUSSIAN PROCESSES 3 



distribution as (\t\ 2H + Z s ( f, . . . , \t\ 2H + Z^) under P, by comparing generating functions: for 
any/3i,...,/3 fc GM, 



logE^exp (E&^J 

= -\t\ 2H -J2^\ 2H + Var 



2 PiCov(B t ,B Si )-Y,Pi\*i\ 2H + Var 



E# s « 



Eft 0* 

i 

£ft|t| 2H +E 



2H_| s ._ t |2Hl +Var 



E# 5 « 



E^ } 



+ - Var 
2 



E 



+ -Var 
2 



E^? 

EA(lf H + 4 } ) 



The first claim of the lemma then immediately follows from the Cramer- Wold device. 

Alternatively, one could carefully define a space on which the distribution of Z becomes a 
Gaussian measure and then note that the claim follows from the Cameron-Martin formula; see 
\\.2\ Prop. 2.4.2] and [20J for key ingredients for this approach. 

When the functional F is translation-invariant, we conclude that 



E Q F(Z) = EF(\t\ 2H + Z {t) ) = EF(Z W - V2B t ) = EF(6 t Z) 
and this proves the second claim in the lemma. 



□ 



The next corollary readily implies subadditivity of E [supo<t<r e Zt ~\ a s a function of T, a 
well-known fact that immediately yields the existence of the limit in (1.1). Evidently we must 



work under the usual separability conditions, which ensure that the supremum functional is 
measurable. 

Corollary 1. For any a < b, we have 









sup e Zt 




r r b 




" pb—a 


E 


sup e * 


= E 


, E 


/ e Zt dt 


= E 


/ e Zt dt 




a<t<b 




0<t<b-a 




.J a 




Jo 



Proof. Applying Lemma[T]for t = a to the translation-invariant functionals F given by F(z) = 
sup a<s<b e Zs ~~ Za and F(z) = e Zt ~ Za dt yields the claims. (The second claim is also immediate 
from~Ee Zt = 1.) □ 



Corollary 2. For T > ; we have 
1 



(2.1) 



-E 



sup e Zt 




= f E 


0<t<T 


Jo 



su P-uT<s<(l-u)T e 
/l 1 uT " )T eZsds 



du. 



Proof. Applying Lemma [T] to the translation- invariant functional F given by 



F(z) 



su Pte[o,T] e 



Jq e Zu du 
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yields that, for any T > 0, 



1 



-E 



sup e 

0<t<T 



1 

T 
1 

T 



E 



E 



f T e z °ds 
sup_ t < s < T _ f 



eft 



and the statement of the lemma follows after a change of variable. 



□ 



The left-hand side of the identity (2.1) converges to 7-L a by definition. The next proposition 



shows that the right-hand side of (2.1) converges to our new representation, thereby proving 



Theorem [TJ The proof of the proposition itself is deferred to Appendix [Aj 
Proposition 1. For any u £ (0, 1), we have 



lim E 

T-»oo 



sn P-uT<s<(l-u)T e ' 
I ( -uT )TeZsds 



E 



Zi 



r^ Zt dt 

J — oo 



< OO. 



Moreover, 



H a = lim / E 

T->oo 



sn P-uT<s<(l-u)T e 
I [ -uT )T eZsds 



du = E 



J— OO 



Apart from establishing Theorem [TJ this proposition gives two ways of approximating 7~Lct • 
The speed at which the prelimits tend to T~L a is different for these two representations. For the 
second 'integral' representation, which is the classical representation in view of Corollary [2j 
the speed of convergence to 7i a can be expected to be slow. Indeed, it is known to be of order 
1/Vr in the Brownian motion case (e.g., [J7]). This is in stark contrast with the speed of 
convergence in the first representation (e.g., for u = 1/2), as analyzed in the next section. Our 
study shows that the slow convergence speed in the classical definition is due to values of u 



close to the endpoints of the integration interval [0, 1] in the right-hand side of (2.1). 

It is instructive to compare our new representation of T~L a with the classical representation 
of Pickands' constant through a discussion of variances. Note that Ee Zs = 1, Vare Zs = 
g Var(z s ) _ e -Var(z s )^ gQ fogfc fa e variance blows up as s grows large. As a result, one can expect 
that sup <£<r e Zt has high variance for large T. Moreover, significant contributions to its 
expecation come from values of t close to T. These two observations explain why it is hard to 
reliably estimate Pickands' constant from the classical definition. 

Our new representation does not have these drawbacks. Let us focus on the special case 
a = 2, for which it is known that %2 = l/V^- Writing iV for a standard normal random 
variable, we obtain that 



n 2 = E 



V2tN-t 2 



Lu^ 2tN - t2 dt 



E 



suptgR e 



-(t-N/V2) 2 



1 



1 



-*dt 



Jr e 

2, so we can expect it to 



It follows from this calculation that M/S has zero variance for a 
have very low variance for values of a close to 2. 

We next present an alternative representation for H a in the spirit of Theorem [TJ The proof 
of Corollary [2] shows that for any locally finite measure fj,, 

r ~i 1 r z 

1 _ 7. f „ su P-uT<s<{l-u)T e S 



(2.2) 



-E 







sup e 


= f E 


0<t<T 


Jo 



fi (T) {du), 



where ^S T \du) = fi(Tdu) /T. Of particular interest is the case where \x is the counting measure 
on r/Z. Then fi^ 7 " 1 converges weakly to Leb/77, where Leb stands for Lebesgue measure. In 
view of this observation, the following analog of Proposition [TJ is natural. The proof is given 
in Appendix [Aj 
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Proposition 2. For any 7] > 0, we have 

n a = E 



supt e R e 



This identity is particularly noteworthy since the integral in the denominator in the repre- 
sentation of Theorem [I] can apparently be replaced with an approximating sum. For a = 2, 
this means that for any 77 > 0, 

dy 



E 



,kyr) 2 — k 2 rj 2 



2. 



We have not been able to verify this intriguing equality directly, but numerical experiments 
suggest that this identity indeed holds. 

We conclude this section with two further related results. For 77 > 0, define the 'discretized' 
Pickands constant through 



lim — E 

T->oo T 



sup e 

k£Z:0<kri<T 



The proof of the next proposition requires discrete analogs of Corollary [2] and Proposition [TJ 
with suprema taken over a grid and integrals replaced by sums (for the first equality). The 
proof is omitted since it follows the proofs of these results verbatim. 



Proposition 3. For any ij > 0, we have 



kei 



E 



j —00 



'/ 



The second representation for Ha in this proposition immediately shows that H a = lim^o % 
by the monotone convergence theorem and sample path continuity. 

A different application of Lemma [l] yields further representations for T-L a and T~L a ■ Let Ft be 
the indicator of the event that the supremum of its (sample path) argument occurs at t. Since 
E[F t (Z)F s (Z)} = for all s ^ t, we have 



E 



sup e 

fceZ:0<fo7<T 



V E [e Zfc 



Fe v (Z) 




SUp Z kr , 
k£Z:-i<k<T/r]-e 



where we use Lemma [T] to obtain the last equality. This can be written as 



T 



sup e 

fceZ:0<fcrj<T 



sup Z krl 
-uT<k7]<{l-u)T 



n {T) (du), 



where, as before, fi^ T \du) = fi(Tdu)/T and [i is the counting measure on r/Z. Note the 
similarity with ( |2.2| ). Taking the limit as T — > 00 requires verifications similar to those in the 
proof of Proposition [2j the details are given in Appendix [Aj The resulting representation is a 
two-sided version of the Hiisler-Albin-Choi representation [31 [23] , and appears to be new. 



Proposition 4. For r] > 0, we have 
and therefore 



lim 77 



SUp Z kr , 

fcez 



SUp Z krl = 

k&L 



From the point of view of simulation, one difficulty with this representation is that one would 
have to estimate small probabilities when rj is small. Unless one develops special techniques, 
it would require many simulation replications to reliably estimate these probabilities. As 
discussed below, such a task is computationally extremely intensive. 
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3. An auxiliary bound 

This section presents a simple auxiliary bound which plays a key role in the next section. 
To formulate it, let be the following approximation of Z% on a grid with mesh rj > 0: 

z v = [ Z vWvl for 1 > °> 
* Iz^^/^ otherwise, 

and define similarly in terms of Bf. 

Let J be a fixed compact closed interval, assumed to be fixed throughout this section. We 
write 

A(r?) = sup(Z t - 2%), 5( V ) = V2 S up(B t - B?). 

teJ teJ 

Define Mj = sup ugJ e z " and S 1 } = fje z ^du. Note that 

^< e A(,)M<i e A(.) 

s v j~ sj-r, e ■ 

Given an event E, we have for r > e EA ^\ 
E(Mj/S];E) 

< E(Mj/S v j; Mj/S 1 } > r/r,) + T -¥{E) 

1 f'°° T T 

= - / P(Mj/5] > y/ V )dy + -P(Mj/S] > t/tj) + -P(E) 

1 r°° T T 

< - P(e A(,?) > y)dy + -P(e A W > r) + -F(E) 

< 1 f (- M^V ) iy + rip (- w^) 1 ) + r P(E) , 

where the last inequality uses BorelPs inequality (e.g., [21 Thm. 2.1.1]). We can bound this 
further by bounding EA(jj). After setting 

K [rj) = sup(Var(Zt) - Var(Z^)), 

we obtain that A(t/) < k(t/) + 5(j]). We next want to apply Theorem 1.3.3 of Adler and 
Taylor [2] to bound E,5(r/), but the statement of this theorem contains an unspecified constant. 
Our numerical experiments require that all constants be explicit, and therefore we directly 
work with the bound derived in the proof of this theorem. Choose r = l/(2i] a / 2 ), and set 
Nj = | J\r^ H . The proof of this theorem shows that 



/ o 

ES{V) - V lo^2) E 2 3/2 r-^yiog(2i+iiV2) =: £(„), 



which is readily evaluated numerically. 

As a result, whenever r > e £ W +K ^, we have 



+ r-p(-M±^z£W£) + r P(B) . 

To apply this bound, one needs to select r appropriately. Note that we may let r depend on 
the interval J. 
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4. Estimation 



This section studies the effect of truncation and discretization of Z on 7~L a . The bounds we 
develop are used in the next section, where we perform a simulation study in order to estimate 

In addition to 7i a and T~L a , the following quantities play a key role throughout the remainder 
of this paper: 



U a (T) = E 



sup_ T < t < T e 
f^ T e Zt dt 



H a (T) = E 



^9-T/r)<k<T/r) 



T/r)<k<T/r) 



where it is implicit that t is a continuous-time parameter and k only takes integer values. 
Throughout, we assume that the truncation horizon T > and mesh size r] are fixed. We also 
assume for convenience that T is an integer multiple of r\. 

We now introduce some convenient abbreviations. For fixed < a\ < a,2 < • • • , we write 



Jo = (— a\,a\) and Jj - 
this section, we use a± 



Jl 



[a,j,aj+ij 



j: 



-a~ 



■j+l, 



T. Write Mj 



j] with j > 1. Throughout 

fj f&dt, M] = sup fc:fcr)eJj e z ^, 



and Sj 



vEk:kr,eJ j e Zk «, and set M = sup^M,-, S = Ejez^j, and S" = £ ieZ S?. The 



length of an interval Jj is denoted by \Jj\- 

The first step in our error analysis is a detailed comparison of H a = E(M/ S 71 ) and E(Mo / Sq ) , 
which entails truncation of the horizon over which the supremum and sum are taken. As a 
second step, we compare E(Mo/<Sg) to H a (T) = E(Mq/5q), which entails approximating the 
maximum on a discrete mesh. 



4.1. Truncation. This subsection derives upper and lower bounds on E(M/5 T? ) in terms of 
E(Mo/Sq). For convenience we derive our error bounds for a 3 - = T{\+^)i~ l for j > 1, for some 
7 > 0. Presumably sharper error bounds can be given when the choice of the a 3 - is optimized. 

4.1.1. An upper bound. We derive an upper bound on E(M/S 7? ) in terms of E(Mo/Sq). Since 
S > Sj for any j £ Z, we have 



E(M/S V ) = E 



Mo 
Sv 



;M = M 



^;M 



M,- 



< 



E(M /^) + ^E 



M 



3 . 



S] 



; Ma > 1 



(4.1) 



<E(M /5 ^) + 2^E 



M 7 - / ~ 

— v2 sup B s > min |s| 



Set £j = {y^maxsgj^. B s > minggj^ |s| Q }. To further bound (4.1), we use the bounds devel- 
oped in Section |H Thus, the next step is to bound W(Ej) from above. We write Tj for r used 



in the j-th. term. Using the facts that B has stationary increments and is self-similar, we find 
that by Theorem 2.8 in pQ, 



E ( max5, 
\seJj 



E 



max B< 

0<s<\Jj\ 



\Jj\ a/2 E ( max B, 



0<s<l 



(4.2) 



< 2\Ji\ a/2 E ( max sN 



0<s<l 



I T-\ a / 2 
\ J j\ , 



where N stands for a standard normal random variable. We derive a bound on P(-Ej) in a 
slightly more general form for later use. It follows from Borell's inequality that, for < a < b, 



(4.3) 



v2 max B s > c + a a I < exp 



sE[a,b] 



[c + a a - V2(b 
46° 
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provided c + a a > y/2(b - a) a / 2 . Specialized to P(£j), we obtain that for j > 1, a" > 
V2(a j+1 - a,) a/2 , 

mri?\^ \ "'j V ' 2[ ''' "' 
f(Ej) < exp < 



[< - V2(a j+1 - a^f \ f (af - 7 «/ W 

— } = exp < - 



4(l + 7 ) c 



provided T > r y2 1 / a . 

Thus, the error is upper bounded by exp(— dT a ) for some constant d as T — > oo. As a result, 
the error decreases to zero much faster than any polynomial, unlike the classical representation 
for which the error can be expected to be polynomial as previously discussed. This is one of 
the key advantages of our new representation. 

4.1.2. A lower bound. We derive a lower bound on E(M/S" ? ) in terms of K(Mq/Sq) as follows: 

M S/? 



E(M/S V ) > E 



qV cV _i_ qV ' — Z-^i 3 



> 



1 + e 



-E 



Mo 
3? 



j¥0 



T ^E(M /^) - ) ' I 



+ e 



Set I? = {eSg < X^j^o ^j?}- ^° a PP^y the technique from Section [3J we seek an upper bound 
on P(E). Let < 5 < T, to be determined later. Since Sq > 77, we obtain 



V'#0 J j>l 

for any probability distribution {qj : j 7^ 0}. We find it convenient to take qj = t/j(l + ip)~^^ /2 
for some ip > and j 7^ 0. An upper bound on Sj for j ' > 1 is 

SJ < [aj + x — aj)e s6J j = 7 aje 'e s6J j . 

For j > 1, we therefore have 

P{S] > er, qj ) < F(e V§nuaf ' SJ i Bs > e V e a ? gj / (jdj)) 

= P(v / 2max J B s > af + log [eqqj / (jaj)} ) 



< 



exp 



log [ eWj /(7«i)] + a? - V2r /2 af 2 



4(1 +j) a af 

provided T is large enough so that the expression inside the square is nonnegative. The last 



inequality follows from (4.3). 



4.2. Approximating the supremum on a mesh. We now find upper and lower bounds on 
E(M /Sj) in terms of E(A^/Sj). 
For the upper bound, we note that 

E(M /^) < e e E(M 7^) + E(M /^; A (r?) > e). 

We use the technique from Section [3] to bound E(Mq /Sq] Ao(t/) > e), which requires a bound 
on P(Ao(fy) > e). Writing Ko(rj) = max(rj a ,T a — (T — rj) a ), we use the self-similarity in 
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conjunction with Borell's inequality and (4.2) to deduce that 



'(Ao(??)>e) < 



< 



< 



< 



sup V2(B t - B?) > e - sofa) 

x te(-T,T) 

■P [ y/2 sup rj a/2 B t > e - kq(t}) ] 
V *e(o,i) / 

■P ( sup B t > '-"oW ) 



exp 



1 



provided e > ko( ? ?)- 

A lower bound on E(Mo/Sg) in terms of E(Mq/Sq) follows trivially: 

E(M /^)>E(M 7^). 

4.3. Conclusions. We summarize the bounds we have obtained. For any e > 0, we have 
derived the following upper bound: 



(4.4) n a < e e E{M$/S2) + E 



Mo 

qV 
°0 



;A (r?) > e 



+ 2^E 



sup -B s > min | s 



«e J, 



6 J, 



where the second and third terms are bounded further using Section [3j Note that this requires 
selecting a r for each of the terms; we will come back to this in the next section. 
For any e > 0, we have derived the following lower bound: 



(4.5) 



l + e 



l + e 



-E 



Mo 
an ■ 

°0 



and we again use Section [3j We note that we may choose a different e for the upper bound 
and the lower bound, which we find useful in the next section. 

5. Numerical experiments 

This section consists of two parts. The first part studies Ha(T) for suitable choices of the 
simulation horizon T and the discretization mesh r/, and uses the previous section to estimate 
bounds on T~L a . In the second part of this section, we present a heuristic method for obtaining 
sharper estimates for H a . 

Simulation of fractional Brownian motion is highly nontrivial, but there exists a vast body of 
literature on the topic. The fastest available algorithms simulate the process on an equispaced 
grid, by simulation of the (stationary) increment process, which often called fractional Gaussian 
noise. We use the method of Davies and Harte |14j for simulating n points of a fractional 
Gaussian noise. This method requires that n be a power of two. In this approach, the covariance 
matrix is embedded in a so-called circulant matrix, for which the eigenvalues can easily be 
computed. The algorithm relies on the Fast Fourier Transform (FFT) for maximum efficiency; 
the computational effort is of order n log n for a sample size of length n. For more details on 
simulation of fractional Brownian motion, we refer to [211. 



5.1. Confidence intervals. Our next aim is to give a point estimate for Ha(T) and use the 
upper and lower bounds from the previous section to obtain an interval estimate for Pickands' 
constant H a . 

The truncation and discretization errors both critically depend on a, but we choose T and 
7} to be fixed throughout our experiments in order to use a simulation technique known as 
common random numbers. This means that the same stream of (pseudo)random numbers is 
used for all values of a. By choosing T and rj independent of a, the realizations of fractional 
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Brownian motion in the n-th simulation replication are perfectly dependent for different values 
of a. As a result, our estimate of % a as a function of a is smoothened without any statistical 
sacrifice. 

Since T and rj are fixed, our estimates for T-L a are likely to be far off from T~La(T) for small a. 
In that regime our algorithm becomes unreliable, since the truncation horizon would have to 
grow so large that it requires ever more computing power to produce an estimate. Any method 
that relies on truncating the simulation horizon suffers from this problem, and it seems unlikely 
that truncation can be avoided. There is some understanding of the asymptotic behavior of 
T-L a as a | |32(, [22] so this regime is arguably less interesting from a simulation point of 
view. Since we cannot trust the simulation output for small a, we focus our experiments on 
a > 7/10. 

Somewhat arbitrarily, we chose to calibrate errors using a = 1, so that our estimates of 
7i a (T) are close to 7i a for a > 1. The closer one sets the calibration point to 0, the higher one 
has to choose T (and thus more computing power). We estimate Ha(T) using 1500 simulation 
replications, which takes about three days on a modern computer for each value of a. We 
carry out the simulation for a = 14/20, 15/20, . . . , 40/20, and interpolate linearly between the 
simulated points. A high-performance computing environment is used to run the experiments 
in parallel. 

We choose the parameters so that the simulated error bounds from the previous section 
yield an error of approximately 3% for a = 1. The most crucial parameter in the error analysis 
is e. We note that a different e can be used for the lower and upper bound, and that e may 
depend on a, so we take advantage of this extra flexibility to carefully select e. For the upper 
bound in (4.4) we use e = 0.005 + 0.025 • (2 — a), and for the lower bound in (4.5) we use 
e = (0.005 + 0.025 • (2 - a))/3. We use T = 128 and rj = 1/2 18 . 

We next discuss how we have chosen the other parameters in the error analysis from Sec- 
tion |4j These have been somewhat optimized. Equation (4.4) and (4.5) produce bounds on 
n oi in terms of Ha (T) in view of Section |3j but this requires selecting some r for each term 
for which Section [3] is applied. We use Tj = 1.3 • (1.005) J_1 for the j-term in the infinite sum, 
and r = 1.4 for any of the other terms. We set 7 = 0.025 for the growth rate of a,j, and we use 
tp = 0.3 for the decay rate of qj. For these parameter values, all event-independent terms in 
Section [3] are negligible. Finally, we replace % a (T) in the resulting bounds with its estimate. 

In Figure [l] we plot our estimates of T-L a (T) as a function of a (blue, solid), along with 
their 95% confidence interval (green, dotted) and our bounds for H a (red, dash-dotted). The 
numerical values are given in Appendix [Bj Note that the errors we find for a < 1 are so large 
that our error bounds are essentially useless. We do believe that the simulated values are 
reliable approximations to 7i a , but the bounds from our error analysis are too loose. 

A well-known conjecture states that T-L a = l/r(l/a) [18], but (to our knowledge) it lacks 
any foundation other than that lim ce j,o'Ha = lim^o l/r(l/a) = 0, %\ = l/r(l), and H2 = 
l/r(l/2). A referee communicated to us that this conjecture is due to K. Breitung. Our 
simulation gives strong evidence that this conjecture is not correct: the function \/T{l/a) 
is the magenta, dashed curve in Figure [TJ and we see that the confidence interval and error 
bounds are well above the curve for a in the range 1.6 — 1.8. Note that we cannot exclude that 
this conjecture holds, since our error bounds are based on Monte Carlo experiments. However, 
this formula arguably serves as a reasonable approximation for a > 1. 



5.2. A regression-based approach. In the previous subsection, we approximated T~L a by 
Wa(T). The main contribution to the error is the discretization step, so we now focus on a 
refined approximation based on the behavior of % a (T) as 77 J, 0. 

This approach relies on the rate at which T-Cl(T) converges to H a (T). We state this as a 
conjecture, it is outside the scope of the current paper to (attempt to) prove it. 

Conjecture 1. For fixed T > 0, we have lim^o V~ a/2 [Ha(T) - W2(T)] G (0,oo). We also 
have lim^o r)- a / 2 [H a - Hi] E (0, 00). 
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Figure 1. Point estimates (blue, solid) and interval estimates (green, dotted) 
for Tia(T) as a function of a. Our error analysis shows that \Ha(T) — H a \ is at 
most 0.03 for a > 1 (red, dash-dotted). We also plot l/T{\/a) (magenta, 
dashed) . 

We motivate this conjecture as follows. We focus on the last part of the conjecture for 
brevity. Since H a = E[M/S V ] by Proposition § we obtain that 



r,- a ^[H a - HI] =¥,\r ] - al \e m - m ' - 1) x 
where = log M 11 and m = log M. The right-hand side equals approximately 



E 



77 a l 2 {m — m r , 



This expectation involves a product of two random variables. The random variable M I? /S" ? 
converges almost surely to the finite random variable Mj S as rj \. 0. Although we are not aware 
of any existing results on the behavior of i]~ a ^ 2 (m — ml 1 ) or its expectation, we expect that the 
random variable r^ OL l 2 {m — vnP) converges in distribution. Indeed, this is suggested by prior 
work on related problems, see Asmussen et al. [5] for the case a = 1 and Hiisler, Piterbarg, and 
Seleznjev |24} [31] for general results on interpolation approximations for Gaussian processes 
(which is different but related). The rate of convergence of mP to m (or for finite- horizon 
analogs) seems to be of general fundamental interest, but falls outside the scope of this paper. 
Conjecture [T] implies that for some c = c(T) > 0, for small 77, we have approximately 

Ul(T)=U a (T)-cr } a l 2 . 

This allows us to perform an ordinary linear regression to simultaneously estimate c and % a {T) 
from (noisy) estimates of Ha(T) for different (small) values of 77 and fixed a. One could use the 
same simulated fractional Brownian motion trace for different values of 77, but it is also possible 
to use independent simulation experiments for different values of r\. The latter approach is 
computationally less efficient, but it has the advantage that classical regression theory becomes 
available for constructing confidence intervals of 7-L a (T). 

Even though we do not have a formal justification for this approach, we have carried out 
regressions with the same simulated trace for different values of 77. The results are reported in 
Figure [2} The simulation experiments are exactly the same as those underlying Figure [TJ and 
in particular we have used the same parameter values. The red, dashed curves are estimates 
for 7ia(T) for r/ = 2~ 14 , 2~ 13 , 2~ 12 , 2 . Using the regression approach, we estimate Ha(T) 
for 77 = 2~ 18 and compare it with our simulation estimates for the same value of 77 (blue, solid). 
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Figure 2. Estimation of rla(T) for different values of rj. 

The two resulting curves are indistiguishable in Figure [2J and the difference is of order 1CP 3 . 
We have also plotted our regression-based estimate of H a (T) (green, dash-dotted). 

It is instructive to look at the resulting estimate for T-L\(T), since we know that rl\ = 1. 
Our estimate for T-L\{T) is 0.9962650, which is indeed closer to its true value. As the number 
of simulation replications increases, we expect much more improvement. 
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Appendix A. Proofs 



Proof of Proposition^ Our proof of (A.l) uses several ideas that are similar to those in Sec- 
tions [3] and [4j so our exposition here is concise. We fix some rj for which l/rj is a (large) integer; 
its exact value is irrelevant. Recall that the quantities Jj, Mj, Sj, M 7 - , S 1 - from Section Shave 
been introduced with respect to parameters < a\ < 02 < • • • . Here we use different choices: 
01 = [2 1 /( 2a )] , aj = aj-x + 1 for j > 2. 



Abusing notation slightly, we write M^_ uT ^i_ u ^ = sup sg [ 



-uT,(l-u)T] e " and S[- uTj (i- u )Tl 



f_ u rp U e B ds. Since M[_ u r,(i-«)r] ~~ ► M and S\- u tji- u )t] ~^ S almost surely as T — > 00 for 
u G (0, 1), both claims follow after showing that 



(A.l) 



lim sup sup IE 

A ^°°T>0«G(0,1) 



' M [-uT,{l-u)T\ M[-uT,(l^u)T] ' 
. &[-uT,(l-u)T] a [-uT,(l— u)T] 
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Write kj = Kj(rj) = sup teJ _. [Var(Zf) — Var(Z^)]. First suppose that — uT and (1 — u)T lie 
in {. . . , — a2, — ai, ai, a2, • • •}• On the event {Mr.^ n^jo = Mj} for some j G Z, we have 

Mr 



<[-uT,(l-u)T] 
'[-uT,(l-u)T] 



< 



2^sup. 6J |B s -B?|+ Kj M i 



(A.2) 



Note that this bound remains valid if —uT and (1 — u)T fail to li e in {. . . , — a%, — oi, ai, 02, • • •} 
Since E |>/2sup. eJo |B,-B?f 



lim 

A— s>oo 



< 00 by Borell's inequality, ( A.l ) follows after we establish that 

0. 



To this end, we observe that for j ' > 1 

2-/2sup seJ . \B s -B^\+k 




3 ;Mj > 1 



4v / 2sup seJ . |B s -B?|+2k.; 



P(Mj > 1) 



4V2sup se[(U] |B S -Bj| 



sup 

teJj 



V2B t > a? 



(a| - V2E[su PtgJj 
8(aj + 1)° 

8(«i + l) a j ' 



where C denotes some constant and we have used (4.2) to obtain the last inequality. Note 
that af > V2 for our choice of aj. The resulting expression is summable, which establishes 
the required inequality by the monotone convergence theorem. □ 



Proof of Proposition^ Our starting point is (2.2) and the accompanying remarks. By The- 
orem 1.5.5 in Billingsley [H], it suffices to show that Leb(-E) = 0, where E consists of all 
u G [0, 1] for which 



lim E 

T^oo 



Ml 



[-u T T,(l-u T )T] 





~ M~ 


= E 


.S" 7 . 





-« T T,(1-« T )T] 

fails to hold for some {ut} with ut —> u. With minor modifications to the bound (A.2) since 
we work with S v instead of S, the proof of Proposition [T] shows that 



lim sup sup E 

A ^°°T>o u e(o,i) 



M [-uT,(l-u)T] _ M[-uT,(l-u)T] 
D [-uT,(l-u)T] D [-uT,{l-u)T] 



>A 



□ 



This implies that E C {0, 1}, so its Lebesgue measure is zero. 

Proof of Proposition^ As in the proof of Proposition [2j it suffices to show that, whenever 
ut — y u G (0, 1), 



lim 1 

T->oo 



,k& 



sup Z kr] 
-u T T<kr)<(l-u T )T 







sup Z kr} 







A sandwich argument readily establishes that 

lim sup Zk v = sup Zkrj- 

T-*°° keZ:-u T T<kri<(l— u T )T 

The claim follows since almost sure convergence implies convergence in distribution. 



□ 
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Appendix B. Simulated values 

This appendix lists our estimates for T-L a (T) in tabular form for rj = 1/2 18 and T = 128, 
along with the sample standard deviation. We also list the lower and upper bounds on % a , 
where we note that these are estimated values since they depend on % a (T) . We cannot report 
these bounds for a < 1, since our choice of parameter values causes the methodology to break 
down. Our methods can be applied with different parameter values to obtain bounds in this 
regime, but this requires more computing time and is not pursued in this paper. 



Q 


estimate ri a \l ) 


sample stddev Mq/Sq 


lower bound Ha 


upper bound H a 


0.700 


1.1888337 


0.5998979 


— 


— 


0.750 


1.1543904 


0.5614484 


— 


— 


0.800 


1.1184290 


0.5257466 


— 


— 


0.850 


1.0855732 


0.4919238 


— 


— 


0.900 


1.0539127 


0.4625016 


— 


— 


0.950 


1.0235620 


0.4360272 


— 


— 


1.000 


0.9946978 


0.4116689 


0.9837218 


1.0250320 


1.050 


0.9674279 


0.3892142 


0.9582444 


0.9956451 


1.100 


0.9424383 


0.3665194 


0.9338777 


0.9687150 


1.150 


0.9191131 


0.3442997 


0.9111406 


0.9435593 


1.200 


0.8963231 


0.3239746 


0.8889154 


0.9190136 


1.250 


0.8743162 


0.3048379 


0.8674489 


0.8953298 


1.300 


0.8532731 


0.2864521 


0.8469212 


0.8726894 


1.350 


0.8322652 


0.2698805 


0.8264114 


0.8501401 


1.400 


0.8121016 


0.2540026 


0.8067235 


0.8285072 


1.450 


0.7922732 


0.2390896 


0.7873523 


0.8072685 


1.500 


0.7727308 


0.2248372 


0.7682494 


0.7863726 


1.550 


0.7531251 


0.2112524 


0.7490677 


0.7654634 


1.600 


0.7342039 


0.1970492 


0.7305511 


0.7453000 


1.650 


0.7155531 


0.1821118 


0.7122884 


0.7254599 


1.700 


0.6970209 


0.1665167 


0.6941287 


0.7057883 


1.750 


0.6782065 


0.1503939 


0.6756727 


0.6858794 


1.800 


0.6585134 


0.1339708 


0.6563256 


0.6651316 


1.850 


0.6384329 


0.1156335 


0.6365762 


0.6440437 


1.900 


0.6176244 


0.0953090 


0.6160842 


0.6222740 


1.950 


0.5944161 


0.0698590 


0.5931803 


0.5981428 


1.998 


0.5663460 


0.0146697 


0.5653943 


0.5692133 
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